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ABSTRACT 

One of the longstanding unsolved problems of planet formation is how solid bodies of a few decime- 
ters in size can "stick" to form large planetesimals. This is known as the "meter size barrier" . In recent 
years it has become increasingly clear that some form of "particle trapping" must have played a role 
in overcoming the meter size barrier. Particles can be trapped in long-lived local pressure maxima, 
such as those in anticyclonic vortices, zonal flows or those believed to occur near ice lines or at dead 
zone boundaries. Such pressure traps are the ideal sites for the formation of planetesimals and small 
planetary embryos. Moreover, they likely produce large quantities of such bodies in a small region, 
making it likely that subsequent N-body evolution may lead to even larger planetary embryos. The 
goal of this Letter is to show that this indeed happens, and to study how efficient it is. In particular, 
we wish to find out if rocky/icy bodies as large as 10 can form within 1 Myr, since such bodies 
are the precursors of gas giant planets in the core accretion scenario. 

Subject headings: planets and satellites: formation — protoplanetary disks — planet-disk interactions 
— methods: numerical 



1. INTRODUCTION 

There are a number of major unsolved problems in 
the theory of rocky planet formation. One is the in- 
famo us "meter size barrier" problem (|Blum fc WurmI 
120081 ) : When dust aggregates reach sizes beyond roughly 
1 cm, they partially decouple dynamically from the 
gas and can reach relative v elocities with respect to 
other particles of up to 50 m/s (jWeidenschilling &: Cuzzil 
Il993t ). Any such collision will likely lead to fragmen- 
tation and/or erosion rather than growth. In addition, 
these bodies tend to drift radially inward toward the 
star at a high speed, and thus quickly get lo st to the 
planet formation process (jBrauer et al.l l200"8[ ) . More- 
over, as laboratory experiments and numerical modeling 
ijZsom et al.ll20ld( ) show, already at millimeter sizes the 
dust aggregates stick insufficiently well for coagulation 
to continue. Despite many theoretical studies aimed at 
solving this problem, this issue is still wide open. In 
recent years a new line of thought has emerged which 
invokes various local partic le trapping mechanisms to 
solve it. iCuzzi "eFall (120081) propose that particle con- 
centrations in a turbulent disk occurring naturally at 
small eddy scales can, statistically, sometimes lead to 
self-gravitating "sandpiles" that gradually condense to 
form planetesimals of 1 to 100 km radius. The particles 
must have grown to millimeter sizes before this process 
sets i n, but the mete r size barrier is thus easily over- 
come. iJohansen et al.l ()2007f ) showed that particle trap- 
ping and subsequent gravitational contraction of parti- 
cle clouds can also happen at the scale of the largest 
turbulent eddies, leading to bodies between 100 km and 
1000 km in size. This mechanism can in fact also act on 
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scales much larger than the turbulent eddy scale. For in- 
stance , iBarge fc Sommerial (|1995[ ) and lKlahr &: Hennind 
(|1997[) showed that particles tend to get trapped in an- 
ticy clonic vortices, i f they exist in protoplaneta ry disks. 
IVarniere fc Taggerl ([2006*) and iTergu cm (2008) showed 
with alpha disk models (Shakura fc Su nvaev 1973) that 
pressure enhancements are expected at dead zone bound- 
aries because of the difference in turbulent viscosity. 
Such enh ancements were co n firmed in the MHD simu- 
lations of iKato et all (|2010D . iKretke fc Li3 ()2007D sug- 
gested that similar pressure ri dges form at the snow lin e 
in th e disk (see discussion by iDzvurkevich et all I2010D , 
while IJohansen et al.l (j2009() showed that an "inverse cas- 
cade" of magnetorotationally driven turbulence will lead 
to large scale pressure bumps in disks. 

Another extensively studied problem is the retention 
of protoplanets in the disk for a sufficiently long time 
for the core accretion process. This requires at least 
10 times slower type I migration speed tha n estimated 
from analytical theory (jAlibert et al.l [20051 ). In recent 
years, however, the understa nding of type I migration 
has c hanged considerably fe.g.'Paardekoopcr fc M ellemal 
20061: Fklcv c t al. 2009; Paardekoopcr et al. 2010). As 



Morbidelli et al.l (j2008i) have demonstrated by using a 



proper surface density profile in their hydro simula- 
tions a planet trap appears, which helps forming mas- 
sive rocky /icy bodies and prevents their migration to the 
central star. We will hereafter adopt the (non-standard) 
nomenclature "terrestrial planet" and "planetary core" 
for rocky/icy bodies of mass below and above 10 M^, 
respectively. 

Many o f the above ideas were combined in a single 
model bv iLvra et aD (|2008L 120091 ). The model explores 
what happens at dead zone boundaries, wher e the den- 
sity e nhan cement was predicted by iVarnier e fc TaggeH 
(|2006f) and llnaba fc Bargj (I2006D to be unsta ble to the 
Rossby Wave Instabilitv (iLovelace et al.lll999l ) , which in 
turn leads to the formation of large scale anticyclonic 
vortices. Particles of about cm to meter in size subse- 
quently drift into the vortices and form gravitationally 
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bound clumps of solids, ranging between the masses of 
the Moon and Mars. 

The presen t paper follows up on that work. iLvra et al.l 
(j2008l l2009[ ) explored only the formation mechanism, 
since the hydrodynamical models arc too burdensome 
to follow the evolution for longer than a few thousand 
of years (which is but a small fraction of the lifetime 
of the gas-rich phase of the disk). A large number of 
massive embryos are formed that all initially have semi- 
major axes very close to the radial location of the pres- 
sure bumps where they were formed. The location near 
these "planet factories" are thus quickly overpopulated 
with embryos. Mutual collisions and merging events are 
naturally expected in such scenario. It would be of a 
great interest to determine how the ensemble develops 
over a long time scale. 

In this paper wc follow the N-body evolution of the 
heaviest embryos of the swarm that were produced dur- 
ing the first few hundred orbits in the above model. We 
include the gravitational interaction between the em- 
bryos and the gas in the disk which leads to type I migra- 
tion. W e account fo r it by applying the analytic formulae 
used bv|Lvra et aD ( 2010 ). and developed originally by 
iPaardekooper et al.l ( 2010 ). The model is detailed in the 
next section. 

2. IMPLEMENTATION OF TYPE I MIGRATION 

To integrate numerically the differential equations 
of the gravitational N-body problem we developed a 
Bulirsch-Stocr integrator, which can handle collisions be- 
tween nearby bodies. When the mutual distance between 
the center of mass of two nearby bodies is less than the 
sum of their physical radii, the bodies collide and merge. 
The physical radius of a body is calculated using its mass 
assuming a 2 g/cm'^ bulk density for the whole embryo 
population. The mass and the initial velocity of the 
newly formed body are calculated assuming perfectly in- 
elastic collision using the center of mass approximation. 

In what follows, we describe briefly how the dissipative 
forces for type I migration have been implemented in our 
N-body code knowing the surface density and tempera- 
ture profiles of the disk. 

In isothermal disks the total torque F experienced by 
a body can be written as 



T/Ta = -0.85-a-0.9^, 



(1) 



where a and /3 are the negative of the local surface den- 
sity and temperature (T'(r)) gradients 
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and Fn = {q/h)'^Y,{r)r'^n{r)'^ (see IPaardekooper et al.l 
I2010D . Here q = m/M^, is the body to star mass ra- 
tio, h = H{r)/r is the disk constant aspect ratio {H{r) 
being the disk's vertical thickness), and D, the Keplcrian 
angular velocity. 

Using Equations (1) and (2), the total torque F can be 
easily determined, which enables the calculation of the 
body's radial migration speed as follows: 
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Fig. 1. — The surface density profile plotted with solid red (dark) 
line, the dimensionless torque F/Fg with the dotted green (light) 
line. The zero torque places are practically at the two "peaks", 
and the three "valleys". If F/Fq > outward, if F/Fq < inward 
migration occur. 



body. The migration timescale is t, 
using Equation (3) reads 



migr 
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Note that a negative Tmigr means outward migration. In 
addition to the inward migration, a body also feels the 
strong damping effects of the disk on its orbital eccen- 
tricity and, since the N-body part of our model is fully 
3-D, its orbital inclination. The corresponding damping 



timescales are To, 



(see the exact formulae 



and Tinr, the corre- 



where L = m\/GM^r is the angular momentum of the 



in Ta naka fc WardI (l200l ). 

Knowing the timescales Tmigr, 
spending forces can be implemented in the N-body code. 
In our code we applied the formula of lCresswell &: NelsonI 
(|2008f) . which for the ith body is 

^'^migr 'ccc 'Hnc 

where k is the unit vector in the z— direction and f^.grav 
is the gravitational acceleration induced by all the other 
bodies (the N-body force). 

3. PHYSICAL BACKGROUND OF OUR SIMULATIONS 

As initial conditions of our N-body experiments we 
used the embryo popul ation form ed during the hydrody- 
namical simulations of iLvra et al. (2008, 2009i). In that 
hydro simulation the initial surface density and temper- 
ature profiles followed power law distributions Ti 
with a ^ -0.5, and T ^ with /? = -1. The den- 
sity profile changed considerably during the simulations, 
but due to the local isothermal disk assumption, /3 = — 1 
was kept fixed. For the torque calculations we used the 
azimuthally averaged surface density profile obtained at 
the end of the hydro simulations, see Fig. 1. To keep the 
problem tractable we fix the gas density profile in time as 
it appears in the figure. The dimensionless torque F/Fq 
is also displayed. Notice that it reaches values as high 
as ±10. We caution that Equation (1) was derived for 
smooth disks, and may thus not be unproblematic when 
used for strongly irregular disks like the one in our model. 
But at present this is the only tractable way to treat the 
problem of type I migration of many bodies over a long 
timescale. 
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In Fig. 1 the two density peaks appear at the inner 
and the outer edge of the dead zone, and correspond 
to the large vortices. Since p = Sc^ (p is the verticaly 
integrated pressure and Cg is the local sound speed) , and 
the temperature distribution (T(r)) is fixed, the density 
peaks correspond approximately to pressure maxima. 

The embryos form in narrow regions around the two 
pressure peaks (corresponding to the density peaks of 
Fig. 1), and they are expected to disperse somewhat due 
to mutual gravitational scattering effects. The peaks are 
very nearby to zero-torque radii. They act as planet 
traps ( "planetary convergence zones '0) because inward 
migration occurs for radii larger than the peak's location, 
whereas outward migration occurs for smaller radii than 
that. There are zero torque places in the density valleys 
as well. But, in contrary to the peaks, they have repulsive 
character. 

On the other hand, the embryos form on a rela- 
tively short timescale, during 200 periods (measured 
at Jupiter's orbital distance) consuming practically the 
solid content of the (strongly truncated) disk used in the 
hydro simulations. If we consider a realistic 100-200 AU 
disk, repeated formation of the planetary embryos can 
be expected as long as the large vortices at the dead 
zone edges exist, and the disk contains enough solid ma- 
terial with size between 1 cm - 1 m. These particles are 
very strongly affected by the gas drag, and therefore are 
drifted inward quickly and continuously to the embryo 
forming vortices as long as gas exists in the system. In 
this way the outer part of the disk acts as a large reservoir 
of solid material and due to the strong radial drift the 
large vortices, which are working as "planet factories" 
are replenished continuously with "raw" material. The 
continuous embryo formation in the "planet factory" is 
an important part of our physical assumptions. In our 
N-body runs we do not intend to model this process, for 
feasibility reasons. But to mimic the continuous produc- 
tion of new bodies, we simply inject Mars-mass embryos 
stochastically into the pressure trap region, given a pro- 
duction rate which we take as a parameter of our model. 
The new bodies are assumed to form only at the out- 
ermost pressure trap (the one at the outer edge of the 
dead zone). We do this for two reasons: One is to make 
the proces more "clean" and understandable, focusing on 
a single pressure trap only. The other is that it is rea- 
sonable to assume that the influx of dust from the very 
outer disk regions gets captured by the outer pressure 
trap, thus choking the inner pressure trap off from the 
supply of new material. 

Armed with the above ideas, we performed two long 
term N-body simulations to see how terrestrial planets, 
and even planetary cores can be formed around the zero 
torque places of the disk. 

4. RESULTS AND DISCUSSION 

Since we are interested in the formation of massive 
planets, we considered as initial embryo population of 
our N-body runs only the 10 most massive bodies, 
formed dur i ng on e of the hydrodynamical simulations of 
ILvra et all (|2009l ) in the outermost pressure trap. We 
performed two simulations: In the first one (simulation 



A) a 1/3 Afe bodjEI was assumed to form at irregular 
time intervals according to a Poisson distribution with a 
rate of 2x 10"''' year~^. In the second one (simulation B), 
a 1/6 Mgg body was assumed to form at time intervals 
also following Poisson distribution with a rate 4 x 10^"* 
year~^. The results of simulations A and B are shown 
respectively in Figs. [2] and [3l 

At the beginning of our N-body simulations close ap- 
proaches, scattering events and collisions happened be- 
tween the initial embryos, after which a multiple mean 
motion resonant structure formed. Subsequently, when 
new bodies got inserted into the pressure traps, this res- 
onant structure is perturbed and new scattering events 
and mergers happen, leading to ever larger masses. The 
whole structure of the resonantly interacting planets is 
broadening, because of the ever increasing masses of the 
planets and the increasing number of them. At around 
10^ to 1.7 X 10^ years the planets which have been reso- 
nantly pushed the farthest away from the pressure bump 
reach the edges of the type I migration convergence zone 
(4.4 "^r < 6.6 AU) and start to migrate away from their 
birthplace. For the ones that migrate inward, they get 
again trapped, this time in the innermost migration con- 
vergence zone. 

Through collisions some of the bodies are able to in- 
crease their masses by accreting either the newly formed, 
or the already existing embryos. After 4 x 10^ years-long 
numerical integration, a planetary core of 10 was 
formed both in simulation A and B. Thus, regarding the 
final mass of the planetary core formed, there is no sig- 
nificant difference between the results of simulations A 
and B. Comparing the evolution of the semi-major axes 
in simulations A and B, it can be clearly seen that in 
case A the evolution of the system is "smoother" than 
in case B. The reason for this is that in simulation A a 
smaller number of bodies is involved than in simulation 
B. However, the final mass of the planetary core formed 
in both simulations is almost the same after ~ 4 x 10^ 
years. This means that the most important parameter 
of our simulations is the amount of mass injected during 
a given time. 

We call attention to the fact that, in some aspects, 
the hydro simulations of lMorbidelli et aTl (j2008[) are sim- 
ilar to ours. We therefore briefly compare the methods 

and results of the t wo work s. In that study the 2D 

hydro code FARGO (jMasseti [2000t) has been used, and 
the N-body integrator was implemented in 3D, in which 
the planet's incl i nation has been damped according to 
iTanaka fc WardI lj2004D . During the hydro simulations 
the steep surface density profile halted the migration of 
10 embryos (having each a mass of 1 M^), but the reso- 
nant structure between them survived only temporarily, 
and some of the embryos were able to collide and merge 
without any additional perturbative event. In our case, 
if there would have been no continuous formation of em- 
bryos, the resonant structure would be stable, preventing 
the bodies from further collisions. There are however a 
few major differences between the physical models. The 
static surface density profile used in our study (providing 
a strong constant torque), differs f rom t hat in the hy- 
dro simulations of lMorbidelli et aTl (|2008[ ). They use an 



* The nomenclature convergence zone was suggested by C. Mor- 
dasini, and adopted here. 



^ Correspondin g to the mass of the largest body obtained by 
ILvra et al.l l|2009l 'l. 
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Fig. 2. — Left: Evolution of the semi-major axes during the simulation A. Right: The corresponding overall mass accretion of the embryo 
population. The most massive body is a planetary core with mass 10 Mgj. 





Fig. 3. — The same as in Fig. [2]for siiimlation B. 



evolving 2D disk, in which the back reaction of the bodies 
to the disk influences the local gas density distribution, 
thus the torques are probably non-constant. The torques 
are also weaker, since the surface density distribution is 
apparently much smoother than in our case. If instead 
of a strong static torque, we would use a weaker and 
oscillating one, the protecting resonant structure might 
be destroyed, enabling even a more rapid forma t ion o f 
a planetary corfj^. Compared to iMorbidelli et al.l (|2008f ) 
the novelty of our approach is that we consider the con- 
tinuous formation of embryos assembled at the pressure 
trap, and we study their long-term evolution. We also 
use smaller building blocks, and follow the core formation 
process through a broader mass range of bodies involved 
in it. 

In the core accretion scenario, the formation of gas gi- 
ant planets can only occur if planetary cores are grown 
still in the gas rich phase. As our simulations show, 
the time required for the assembly of a 10 core 
is < 0.5 X 10^ years in the framework of our settings. 
Since the lifetime of the gas disk is expected to be sa 5 
Myr, there is significant time left for accretion of nebu- 
lar gas, co mpleting the formation of a gas giant planet 
(jAlibert et al. 2005). 

® The torque oscillation would be due to the fact that the trap- 
ping zone i n fact consists of a couple of large vortices arranged in 
a circle (see lLvra et al.|[2008l ). which means that, depending on the 
azimuthal location of the planet compared to the locations of those 
vortices, the planet may experience different torques. 



5. SUMMARY AND OPEN QUESTIONS 

In this Letter we investigated the possibility to form 
large bodies in protoplanetary disks at such places where 
the torques responsible for type I migration vanish. We 
assumed that at the edges of the dead zone of the disk 
large vortices develop that can collect the (from centime- 
ter to meter sized) solid content of the disk. Through the 
self gravity of these overdense regions of solids, relatively 
large embry os arc formed, with masses up to 1/3 M0 
(jLvra et al.,.2008, , 2009.) . 

Due to the particular surface density distribution ob- 
tained from the above mentioned hydro simulations the 
embryos stay trapped close to their birthplaces because 
the location of the type I migration convergence point 
is very close to the location of the pressure trap. The 
embryos capture each other into mean motion reso- 
nances, entering into a very robust protective configu- 
ration against further collisions. This resonant structure 
is, however, perturbed from time to time when a new 
massive embryo forms in one of the giant vortices. Dur- 
ing these perturbative events, the embryos can collide 
and, by merging events, form planetary cores as massive 
as 10 M^. The whole process takes < 0.5 Myr. 

We stress that the proposed scenario should also work 
in disk models where the pressure trap em erges by other 
mechanics, such as at the ic eline. as in [ Kretke fc LinI 
(|2007t ): or in zonal flows, as in iJohansen et aLf i 2009f ) . 

Besides cores of giant planets, other massive terrestrial 
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planets can also be formed, which by scattering at later 
epochs may lead to the formation of a complete plan- 
etary system. The escape of the bodies from the zero 
torque places can also be caused by the slowly forming 
giant planet, which by opening a gap will change the sur- 
face density profile. On the flip side, however, this sce- 
nario would predict that, if the trap is located beyond 
the ice line, all rocky planets of the resulting planetary 
system should be ice- or ocean planets, which for our 
solar system is clearly not the case. Further research is 
thus required to investigate these issues. 

It is clear, however, that our present model is still very 
simplified, and further study is required to verify our pro- 
posed scenario for t he formation of rocky/ icy planetary 
cores. For instance iMorbidelli et al.l (|2008[ ) have shown 
that, due to the large vortices at the planet trap, the 
semi-major axes of the embryos oscillate. This means 
that the position of the zero torque point of any given 
planet is not constant in time. Thus the dynamics of the 
whole embryo population is perturbed, which may result 
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in more effective collisions of bodies at the convergence 
zone, shortening the time to form a planetary core. On 
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posed to be formed, thus perhaps quenching the forma- 
tion of new embryos. Also, the planets that are pushed 
away from the convergence zone to larger radial distances 
may form a barrier to dust drifting inward from the outer 
parts of the disk. All of these questions require fur- 
ther study, and much more detailed modeling. But con- 
sidering that the combined dust- and planet-migration 
trap scenario we propose here has the potential to solve 
both the meter-size barrier problem and the time scale 
problem of oligarchic growth, it seems worthwhile to in- 
vest substantial effort in studying planet formation along 
these lines. 
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